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Several simulations of turbulence in the Large Plasma Device (LAPD) [W. Gekelman 
et al, Rev. Sci. Inst. 62, 2875 (1991)] are energetically analyzed and compared 
with each other and with the experiment. The simulations use the same model, 
but different axial boundary conditions. They employ either periodic, zero-value, 
zero-derivative, or sheath axial boundaries. The linear stability physics is different 
between the scenarios because the various boundary conditions allow the drift wave 
instability to access different axial structures, and the sheath boundary simulation 
contains a conducting wall mode instability which is just as unstable as the drift 
waves. Nevertheless, the turbulence in all the simulations is relatively similar because 
it is primarily driven by a robust nonlinear instability that is the same for all cases. 
The nonlinear instability preferentially drives kn = potential energy fluctuations, 
which then three-wave couple to k\\ 7^ potential energy fluctuations in order to 
access the adiabatic response to transfer their energy to kinetic energy fluctuations. 
The turbulence self-organizes to drive this nonlinear instability, which destroys the 
linear eigenmode structures, making the linear instabilities ineffective. 
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I. INTRODUCTION 



Hydrodynamic turbulence often occurs in the absence of linear instability, e.g. turbulence 
in pipe flow (Pouseille flow}^. Although many robust linear instabilities exist in magnetized 
plasmas, nonlinear instability can arise as has been shown in many turbulence simulations. 
Sometimes the instabilities were found to be of the subcritical2~- or supercritical^ variety, 
while at other times the nonlinear instability simply overpowered a particular linear instabil- 
ity to drive the turbulence^— . In the scenario applicable to turbulence in the Large Plasma 
Device (LAPD)— - in which the magnetic fields lines are straight and without shear - several 
studies showed the turbulence to be driven by a nonlinear instability of this last type, where 
the nonlinear instability imposed itself over a linear drift-wave instability^ 1 ^^ 1 ^. One of 
these^ used LAPD experimental parameters and profiles in the simulation and demonstrated 
that the nonlinear instability was necessary to drive turbulence with characteristics similar 
to that of the experiment. Now, these studies^^^ 1 ^ all ascertained that the mechanism 
driving the nonlinear instability relies upon axial wavenumber transfers between ku = and 
k\\ 7^ structures. The reason is that the turbulence self-organizes to preferentially drive 
ku = density and temperature fluctuations, taking energy from the density and tempera- 
ture equilibrium gradients. But in order to access the adiabatic response, which transfers 
energy into the dynamically critical E x B flows, the ku = fluctuations must transfer their 
energy through nonlinear three-wave decay into ku ^ fluctuations. 

All of these straight magnetic field simulations, however, employed periodic boundary 
conditions in the axial (field-aligned) direction. One can justifiably question the use of 
periodic axial boundary conditions to model LAPD and wonder whether the nonlinear in- 
stability is an artifact of this choice. The use of this boundary condition, for example, 
prevents the fastest growing linear drift-waves in LAPD from being captured in the simu- 
lation: with periodic boundary conditions the longest wavelength mode has its wavelength 
equal to the length of the device whereas the fundamental mode (half-wavelength equal to 
the device length) has a higher growth rate. This may be the reason why the simulations 
produce k\\ = structures. Would these disappear if the simulations are allowed access 
to longer wavelength structures through different boundary conditions? Furthermore, the 
actual LAPD axial boundary includes conducting structures; it is well known that sheaths 
on metal walls can drive linear instabilities like the conducting wall mode^. Clearly, pe- 
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riodic boundary simulations miss this wall physics. Axial boundary conditions do in fact 
have a significant affect on the linear stability properties of the system, which is shown 
in detail in Section IHII And since the nonlinear instability relies upon very-long-parallel- 
wavelength modes to operate, it is reasonable to speculate that the bondaries also might 
have a significant affect on the nonlinear instability properties. 

The real axial boundaries in LAPD are complicated. One side of the device contains a 
hot, emissive cathode behind a mesh anode, and in front of that sits a biased limiter with 
radius slightly less than the cathode radius^. The other side contains a floating mesh plate, 
which is likely shielded by a layer of neutral gas (the plasma may be detached from the 
plate on this end). Not only is it difficult to determine what the actual axial boundary 
conditions are, it is also difficult to develop and implement models for the boundaries. 
Thus, this paper takes a simpler approach of exploring the affect of non-periodic axial 
boundary conditions on the nonlinear turbulent dynamics using various idealized boundaries, 
leaving the calculation and implementation of physically realistic boundary conditions to 
future work. The different axial boundary conditions used here include zero-value, zero- 
derivative, and perfectly conducting metal plates, which can be modeled with a Bohm 
sheath condition. The main finding is that the nonlinear instability is robust to changes 
in boundary conditions: while the linear stability properties are modified significantly with 
different boundary conditions, the nonlinear instability still dominates the turbulent drive 
in all cases. In fact, the qualitative properties of the turbulence and the turbulent dynamics 
are similar between the simulations. Quantitatively, there are some differences between the 
simulations such as varying fluctuation levels and varying degrees to which the nonlinear 
instability dominates the linear ones. 

The paper is organized as follows: Section HT1 presents the model and boundary conditions 
used in the simulations, while Sec. IHII goes over the origins and properties of the linear 
instabilities in the different simulations. Section IIVI develops the energetics equations that 
are used in Sec. |V] to reveal the details of the nonlinear instability in the simulations. 
Finally, Sec. I VI I explores the effect of the nonlinear instability on the mode structure of the 
turbulence. 
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II. THE SIMULATION MODEL 



A Braginskii-based fluid model^ 2 - is used to simulate global drift wave turbulence in LAPD 
using the BOUT++ code 23 . The evolved variables in the model are the plasma density, N, 
the electron fluid parallel velocity v\\ e , the potential vorticity w = Vj_ • (N Vj_(fi), and the 
electron temperature T e . The ions are assumed cold in the model (Tj = 0), and sound wave 
effects are neglected = 0). Details of the simulation code, derivations of the model, grid 
convergence studies, and analyses of simplified models may be found in previously published 
LAPD simulation studies^^— . 

The equations are developed with Bohm normalizations: lengths are normalized to the 
ion sound gyroradius, times to the ion cyclotron time, velocities to the sound speed, densities 
to the equilibrium peak density, and electron temperatures and potentials to the equilibrium 
peak electron temperature. These normalizations are constants (not functions of radius) and 
are calculated from these reference values: the magnetic field is 1 kG, the ion unit mass is 
4, the peak density is 2.86 x 10 12 cm -3 , and the peak electron temperature is 6 eV. The 
equations are: 



d t N = -v E • ViV - AWiRie + MjvViiV + S N + {0, N}, (1) 

d t v lle = -— ^VyiV - 1.71^V||T e + ^V|,0 - u e v ]{e + {<f>,v ]le }, (2) 
m e Ao m e m e 

d t zu = -N V\\V\\ e - v in w + [i^sl\w + {<p, w}, (3) 

2 2 
d t T e = -ve • VT e0 - 1.71-TeoV||f|| e + ^^K|| e V|T e 

e T e + fi T V 2 ± T e + S T + {<P,T e }. (4) 
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In these equations, /^t, and fi^ are artificial diffusion and viscosity coefficients used 
for subgrid dissipation. They are large enough to allow saturation and grid convergence 2 ^., 
but small enough to allow for turbulence to develop. In the simulations, they are all given 
the same value of 1.25 x 10 -3 in Bohm- normalized units. This is the only free parameter 
in the simulations. All other parameters such as the electron collisionality z/ e , ion- neutral 
collisionality i/j n , parallel electron thermal conductivity K\\ e , and mass ratio ^ are calculated 
from experimental quantities. There are two sources of free energy: the density gradient due 
to the equilibrium density profile Nq, and the equilibrium electron temperature gradient in 



T e o, both of which are taken from experimental fits. Nq and T £ q are functions of only the 
radial cylindrical coordinate r, and they are shown in Fig. [TJ The mean potential profile 
4>o is set to zero in the model, and terms involving O are n °t included in Eqs. dH The 
justification for this is that biasable azimuthal limiters in LAPD allow for the mean E x B 
flow and flow shear to be varied with high precision, even allowing the flow to be nulled 
out>2i. The simulations in this paper use the N Q , T e0 , and 0o profiles from the nulled out 
flow experiment, justifying setting the mean potential profile to zero in the simulations. 

Simulations also use density and temperature sources (S n and St) in order to keep the 
equilibrium profiles from relaxing away from their experimental shapes. These sources sup- 
press the azimuthal averages (m — component of the density and temperature fluctuations) 
at each time step. The azimuthal average of the potential <fi is allowed to evolve in the simula- 
tion, allowing zonal flows to form, although they are relatively unimportant to the turbulent 
dynamics^. 

The terms in Poisson brackets are the E x B advective nonlinearities, which are the only 
nonlinearities used in the simulations. The numerical simulations use finite differences in all 
three dimensions and use cylindrical annular geometry (12 < r < 40 cm). The radial extent 
used in the simulation encompasses the region where fluctuations are above a few percent 
in the experiment. Therefore, the radial boundaries are fixed to zero value. 

This study analyzes five turbulent simulations which will be referred to as (1) the periodic 
simulation, (2) the n = suppressed simulation (n is the axial wavenumber), (3) the sheath 
simulation, (4) the Dirichlet simulation, and (5) the Neumann simulation. The periodic 
and n = suppressed simulations were also analyzed in a previous paper 1 -. Both of these 
simulations enforce periodic boundary conditions in the axial (z) direction. The n = 
suppressed simulation adds an artificial sink-like contribution to Eqs. [T0 which removes the 
axial average (k\\ = contribution, where n = k\\l\\/2ii) of the fluctuations at each time step. 
This n = suppression eliminates nonlinear instability drive and allows the linear instability 
to take over the turbulent drive^. The n = suppressed simulation, therefore, serves as a 
contrast to the periodic simulation in which a nonlinear instability drives the turbulence. 

The sheath simulation, as its name implies, uses sheath boundary conditions on the 
axial machine ends. Specifically, the sheath boundary condition for the parallel current is a 
linearized Bohm condition: 
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J\\ = N ((f> + logy/iwmjmi T e ) 



(5) 



where J\\ = —N v\\ e . The axial boundary for is set using this relation along with 
Ohm's law: ^-V||0 = v e v\\ e . The axial boundaries for the density and temperature fields are 
implemented with zero-second-derivative boundary conditions. This is somewhat arbitrary, 
and it is noted that stringent analytical and numerical calculations have recently been made 
for such fields in the magnetic pre-sheath region 28 , but those have not been applied in 
this simulation. The fourth (Dirichlet) simulation uses fixed zero-value axial boundary 
conditions, while the fifth (Neumann) employes zero-first-derivative conditions to all fields. 

As a first comparison, Fig. [2] shows a few statistical characteristics of the density fluctu- 
ations for each of the five simulations along with the corresponding characteristics from the 
experiment. The periodic simulation clearly has the most similar characteristics to those of 
the experiment while the n = suppressed simulation is most dissimilar. The fluctuations of 
the sheath, Dirichlet, and Neumann simulations have similar statistical properties as those 
from the experiment, however, the amplitude of the fluctuations is a bit smaller than the 
experimental fluctuations in general. 
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FIG. 1. The profiles of density Nq and electron temperature T e o used in the simulations normalized 
to their peak values of 2.86 x 10 12 cm -3 and 6 eV, respectively. 
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FIG. 2. A comparison of the statistical turbulent properties of the density field. The figures show 
a) the frequency spectrum, b) the pdf and c) the RMS level of the fluctuations as a function of 
radius. The data are calculated from the experiment (black), periodic simulation (red), the sheath 
simulation (blue), the n = suppressed simulation (green), the Dirichlet simulation (magenta), 
and the Neumann simulation (yellow). This color scheme is consistent throughout the paper. 
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III. LINEAR INSTABILITIES 
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FIG. 3. The linear growth rates of the four simulations (the n = suppressed simulation has the 
same linear growth rate as the periodic simulation) along with the growth rate of the conducting 
wall mode (CWM) as a function of azimuthal wavenumber m. 

The model described by the equations of Section [III contains a few linear instabilities that 
can all act at the simulated scales. Two of these instabilities are of the electrostatic drift 
wave type - one is driven by the density gradient and the other by the electron temperature 
gradient. Both of these instabilities supply energy to the electrostatic potential through 
parallel compression, called the adiabatic response, and are made unstable by the electron- 
ion collisional dissipation. These instabilities act under all choices of parallel boundary 
conditions. 

The other instability is called the conducting wall mode (CWM) since it is driven by 
the conducting wall sheaths on the parallel boundaries. Various terms in Eqs. [TH] can be 
eliminated to isolate the conducting wall mode. The reduced set of linearized equations is: 



S 



d t V\\ e = — V||0 - V e V\\ e , (6) 
Tfl e 

d t w = -iV V||U[| e - u in w + /i^Vitu, (7) 

d t T e = -v E • VT e0 + 7^r^||eVjT e - 2 —v e T e + /i T ViT e , (8) 

along with the axial boundary condition given in Eq. The free energy source for this 
instability is the electron temperature gradient, which is also the free energy for the thermally 
driven drift waves. However, the adiabatic response is replaced here with a coupling to the 
potential through the axial boundary condition. 

For the experimental parameters and profiles used in the turbulent simulations, the linear 
growth rates for the drift waves and the CWM are comparable. The growth rates of the 
fastest growing linear eigenmode at a given azimuthal wavenumber are shown in Fig. El The 
periodic, sheath, Dirichlet, and Neumann growth rate curves are found by simulating the 
linearized versions of Eqs. [T]|4]with their respective axial boundary conditions in BOUT++. 
Therefore, both the density- and temperature-driven drift wave contributions are present 
for these curves. The CWM curve is obtained by simulating Eqs. [6j|8] with the sheath axial 
boundary condition of Eq. [5j so there is no drift wave contribution to this curve due to the 
absence of the adiabatic response in these equations. 

The linear growth rates for the Dirichlet and Neumann simulations are markedly different 
from those of the periodic simulation because the zero-value and zero- derivative boundary 
conditions allow for more freedom of the axial structure. In other words, they allow for the 
axial wavenumber n = k\\l\\/2ir to take on non-integer values (about 1/2 in this case) that 
are more unstable than the n — 1 eigenmodes that are enforced by the periodic boundaries. 
The sheath boundary condition has this affect as well, but more importantly it affects the 
linear stability properties of the system at low azimuthal wavenumber m due to the presence 
of the CWM contribution. The CWM is more unstable than the drift waves for m J$ 30 
but less unstable for m <J 30, and the eigenmode with the highest growth rate has m = 20 
for the sheath simulation as opposed to m = 60 for the pure drift wave simulations. The 
linear sheath eigenmodes also reflect which of the linear instabilities is active at which 
wavenumber. In other words, the sheath eigenmodes have CWM character at m J$ 30 and 
drift wave character at m ^ 40. This manifests itself as differences in axial and radial 
structure as well as in phase relations between the different scalar fields. However, the 
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linear differences are only significant in the end if they affect the turbulence dynamics. And 
following a previous paper—, an effective way to study the turbulence dynamics is with an 
energy dynamics analysis of the turbulent simulations. 



IV. ENERGETICS FORMALISM 

In order to perform an energy dynamics analysis on the simulations, expressions for the 
energy and energy evolution must be derived from Eqs. [10 To start, an expression for the 
normalized energy of the wave fluctuations in the model is defined as: 



2jy 



Po ( (A/Ao) 2 + \{T e /T eQ A + N f^vl + (V ± 0) 2 



dV, 



(9) 



where P = N T e0 is the equilibrium pressure. The ^P (N/N ) 2 term is the potential 
energy due to density fluctuations, jPo(T e /T e o) 2 is the electron temperature fluctuation 
potential energy, \No^v 2 e is the parallel electron kinetic energy, and |iVo(V_i_0) 2 is the 
E x B perpendicular kinetic energy. 

A more detailed look at the energetic processes comes from a spectral energy analysis. To 
do this, each fluid field (N, T e , v\\ e , (f>) at a given time is Fourier decomposed as F(r,9,z) = 
Sfc f^(r)e l(jne+kzZ \ where the subscript k represents the spectral wavenumbers, (to, n), and 
both positive and negative wavenumbers are included in the sums, to is the azimuthal 
wavenumber while n is the axial integer wavenumber. Note that the radial direction is not 
spectrally decomposed because it's not essential here. With this, the energy of each Fourier 
k = (to, n) mode is 



2T, 



(0 



rrii 



dr 



(10) 



where the brackets () represent the radial integral: J^ b rdr. The energy evolution for each 
Fourier mode of each field has the form: 



= Q 3 {k) + C s {%) + D 3 {k) + £ k'). 



(11) 



The index j stands for each field, (n,t,v,(f>), and the sum over j gives the total energy 
evolution. Note that with the conventions used, the symbol n denotes both the axial mode 
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number as well as the Fourier coefficient of the density fluctuation. The differences should 
be clear in context. The derivation of Eq. [11] is given in the previous work 18 . Tj(k,k') 
is the nonlinear energy transfer function that comes from the advective nonlinearities. It 
describes the nonlinear energy transfer rate of fluctuations with k' = (m', n') to fluctuations 
with k = (m,n). For example, a positive value of T t (k,k') indicates that temperature 
fluctuations at wavenumber k gain energy from temperature fluctuations at wavenumber k', 
with the process mediated by flow fluctuations at wavenumber k — k'. 

The linear terms are broken up into three contributions in Eq. [TTJ Dj(k) represents 
energy dissipation due to collisions, artificial diffusion and viscosity, and the density and 
temperature sources. Each contribution to Dj{k) is negative. Cj(k) contains the linear 
terms dubbed "transfer channels"—. They are: 



C n {k) = Re{{-ik z T e0 vpil)} (12) 

C v (k) = Re {(-ik z T e0 n^l + ik z N^t - l.mk x N tp)§ } (13) 

CV(fc) = Re{{ik z N Q v n 4>l)} (14) 

C t (k) = Re{(-1.7Uk z N QVf :t*:)} (15) 



First, note that the real part operators are written explicitly in these expressions since 
the imaginary part of these expressions would cancel with the imaginary part of the corre- 
sponding terms with — k. Second, notice that C n {k) +C v (k) + C t f>{k) +C t (k) = 0. This is the 
reason why these terms are called transfer channels. They represent the transfer between 
the different types of energy of the different fields (N, 0, T e v\\ e ), but taken together, they 
do not create or dissipate total energy from the system. The only energy transfer between 
different fields in this system is through the parallel electron velocity (parallel current) dy- 
namics. There is no direct transfer between the state variables N, <f), and T e . Altogether, 
the coupling through the parallel current is called the adiabatic response. It is an essential 
part of both the linear and nonlinear drift wave mechanisms^ 1 ^. The adiabatic response 
moves energy from the pressure fluctuations to the perpendicular flow through the parallel 
current. 

Finally, the Qj(k) terms represent the fluctuation energy sources. They are: 
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Qn(k) 


= Re 


Q v (k) 


= 


qS) 


= 


Qt{k) 


= Re 



' n ' Tea -d r N^l)\ (10) 



N r 



(17) 
(18) 

iW****)} (19) 

Q n (k) is the energy extraction from the equilibrium density profile into the density fluc- 
tuations. This term may have either sign depending on the phase relation between 0^ and 
n^, so it can in fact dissipate fluctuation potential energy from the system as well as create 
it at each k. Qt{k) is completely analogous to Q n (k) but for the temperature rather than the 
density Q v {k) and Q<f>(k) are zero because the parallel and perpendicular flow fluctuations 
obtain energy only through the adiabatic response, not directly through the free energy in 
the equilibrium gradients. 



V. ENERGY DYNAMICS RESULTS 

The diagrams in Fig. H] summarize the flow of energy for the periodic and sheath simula- 
tions. Each of the functions, such as Q n (m, n), is a function of m and n, making visualization 
of all of these functions difficult. So the terms in the diagrams are summed over m. Addi- 
tionally, all of the n^O terms are summed over as well. The n = contribution is separated 
from the other n components because the n = 0f>n/0 dynamic is the primary factor that 
determines whether the linear instability or the nonlinear instability dominates the energy 
drive 18 . 

In these diagrams, the source of energy into the fluctuations is free energy in the equilib- 
rium gradients, ViV and VT e0 . The arrows labeled Q n and Q t represent energy injection 
from the equilibrium gradients into the fluctuations, n(k) and t(k). The four Q arrows con- 
tain values that sum to 100 (by choice of normalization). Since the Q pathways are the only 
pathways that deposit net energy into the fluctuations, the numbers in all arrows represent 
a percentage of the total energy injected into the system. Now, a majority of the energy 
deposited into the fluctuations (71% for the periodic simulation and 56% for the sheath 
simulation) is from the density gradient into the n = density fluctuations. This is not a 
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FIG. 4. Summary of the energy dynamics for the a) periodic case and b) sheath boundary case. 
Each arrow contains the sum over all m. The O's in the parentheses represent the n = modes 
while the !0's represent a sum over the n modes for The values in the arrows represent the 

percentage of total energy that goes through the channel represented by the arrow. The size of 
the arrows gives a rougher but more visual indication of the amount of energy going through each 
channel. 

path allowed by the linear drift-wave instabilities in the system since they can only deposit 
energy into n ^ fluctuations. In fact, in this turbulent state, more energy is transfered 
by nonlinear three-wave coupling into n ^ fluctuations than by direct injection from the 
equilibrium gradients. The three-wave coupling is represented by the T n and T t arrows. The 
direction of these arrows is from ^ = 0411/0, which is opposite to that expected from 
the common cascading type turbulent paradigm where the linear instability dominates the 
turbulent injection dynamics and three-wave processes transfer energy to waves that are 
linearly stable. 

The reason why all of the non-dissipated energy that is injected into n = density and 
temperature structures goes into n^O density and temperature potential energy structures 
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rather than into n = kinetic energy structures is that potential to kinetic energy transfer 
can only work through the adiabatic response, which requires Actually, in the sheath 

simulation, potential energy can transfer to kinetic energy through the axial boundaries, but 
this still requires ti^O and it works only through the temperature fluctuations, which are 
less important than the density fluctuations in the simulations. Note, in fact, that boundary 
contributions aren't included in the energy dynamics calculations as they are insignificant. 
So the main transfer channel from potential to kinetic energy, shown by the C n , Ct, and 
arrows, is the adiabatic response. The adiabatic response takes energy from the n^O density 
and temperature fluctuations and transfers it into the parallel velocity fluctuations (C n and 
C t , respectively). It then transfers some of this energy into n / potential fluctuations 
(C^,), while much of the fluctuation energy is ohmically dissipated (D v ). 

The final step of the energy dynamics process is a three-wave axial wavenumber transfer 
from potential fluctuations at n ^ to potential fluctuations at n — 0, which are of course 
neccessary in making the Q„(0) and Qt(0) terms finite. Meanwhile, dissipation acts on all 
fluctuations, which is quantified by the D arrows throughout. A couple of interesting features 
are evident from the diagrams in Fig. HI First, the turbulent dynamics in both simulations 
are dominated by the nonlinear instability process described above and in Friedman et 
al.— rather than the paradigmatic process of linear instability energy injection followed by 
nonlinear cascading. And second, the periodic and sheath simulations have qualitatively 
similar dynamics despite the fact that the linear stability properties of the two cases are 
qualitatively different. This is also true of the Dirichlet and Neumann simulations, although 
their diagrams are not shown in Fig. HI This speaks to the robustness of the nonlinear 
instability. 

A more compact way to see the similarity of the instability process in all of the simulations 
is shown in Fig. This figure shows the turbulent growth rate of the five simulations. The 
turbulent growth rate is defined as the net energy injection into the fluctuations minus the 
dissipation out of them, all divided by the total energy. The conservative transfers (C's and 
T's) are of course not part of the growth rates. Formally, 



where the index j represents the different fields. Since the growth rates sum over the 




(20) 
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fields, it's not as difficult to view them in their full wavenumber space (both in m and 
n). However, almost all of the energy in the simulations is contained in n = and n — 1 
fluctuations, so n > 1 fluctuations are not shown in the figure. It should be noted, however, 
that Fourier decomposing the non-periodic simulations (sheath, Dirichlet, and Neumann) 
in the axial direction is less ideal than doing so in the periodic simulations (periodic and 
n = suppressed). Fourier modes are not as natural a basis in the non-periodic simulations, 
and the n = 1 Fourier mode does not perfectly represent the linear eigenmode structure 
as it does for the periodic simulations. A more detailed discussion of this point is left to 
the Appendix, but in summary, the n = 1 Fourier mode does capture enough of the linear 
eigenmode structures of the non-periodic simulations to make the Fourier decomposition 
useful for this simulation. 

Figure ^illustrates the true dominance of the nonlinear instability in the generally positive 
n = energy growth rate and generally negative n = 1 growth rate for all of the simulations 
except the n = suppressed simulation. Take the periodic simulation as the clearest example 
of this point since there is no ambiguity in the Fourier transform for this simulation. Fig. [3] 
shows the linear growth rate as a function of m only, but implicitly, n = 1 for the periodic 
curve. The reason is that the linear eigenmodes of the periodic simulation are Fourier 
modes and all eigenmodes with n > 1 have smaller 7zi n (m) than the n = 1 eigenmodes 
and all n = flute eigenmodes have large negative 7z« n ( m ) because the linear instability 
mechanism doesn't work when n = 0. On the other hand, the turbulent growth rate curves 
in Fig. [5] for the periodic simulation have a very different nature than the corresponding 
linear growth rate curves. For low m, the n = turbulent growth rate is positive, while the 
n = 1 turbulent growth rate is negative for all m. This reversal in sign is indicative of the 
nonlinear instability. The sheath, Dirichlet, and Neumann simulations are not as easy to 
analyze because when their most unstable linear structures are axially Fourier decomposed, 
they contain pieces of all n components (their eigenmodes are not Fourier modes). However, 
it's clear that the turbulent growth rates are quite similar between these simulations and 
to the periodic one, indicating that the same nonlinear instability mechanism is acting in 
each of these cases as well. The similarity in the turbulent structures (see Fig. [2]) also 
points to this conclusion. This is in stark contrast to the n = suppressed simulation, 
which has an n = 1 growth rate similar to the linear growth rate (Fig. [3]) and a negative 
n = growth rate. Note that even though the n = fluctuation components are removed 
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from the n = suppressed simulation at each time step, n = fluctuations are nonlinearly 
excited (by three-wave transfer) by the n ^ fluctuations, and therefore, they do have 
small but finite amplitude prior to their removal, which can be used to calculate the growth 
rate of these modes. Furthermore, the turbulent growth rate of the n = 1 component of 
the n = suppressed simulation is slightly less than the linear growth rate due to the fact 
that eigenmodes other than the fastest growing ones are nonlinearly excited in the turbulent 
simulation, thus damping the growth rate. But this n = suppressed growth rate picture 
is just that of the turbulence paradigm of linear instability with cascading dynamics, which 
is significantly different than the nonlinear instability picture of the other simulations. 

0.004 | 1 1 1 1 1 
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FIG. 5. The turbulent growth rates with 7 defined in Eq. [20j The n = (solid) and n = 1 
(dashed) growth rates are displayed as a function of m. The different colors represent different 
simulations, consistent with the scheme used in Fig. [2] 

VI. LINEAR VS NONLINEAR STRUCTURE CORRELATION 

Now it may be the case that in simulations dominated by a linear instability, the fastest 
growing linear eigenmode dominates the system, nonlinearly transfering some energy to more 
weakly unstable or even stable eigenmodes. In this case, a large portion of the energy may 
remain in the fastest growing linear eigenmode^. In the case where a nonlinear instability is 
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dominant, the linear eigenmode should have little bearing on the structure of the turbulence 
and therefore little energy should be contained in this eigenmode. Therefore, a gauge of 
whether a linear or nonlinear instability dominates a system is the fraction of energy in a 
turbulent system that is contained in the fastest growing linear eigenmode. This may be 
calculated by projecting the fastest growing eigenmode onto the turbulent state. 

Formally, in the model considered in this study, the turbulent state is fully described by 
four independent fields, which can be appended into a single vector of the spatio-temporal 
field functions: fturbif^t) = {N(r,t),T e (r,t),'V±4>(r,t),v\\ e (r,t)}. This vector may be de- 
composed in a complete basis: 



where ipi t m{ r , z ) are time-independent spatial complex basis functions of the form ip ijm (r, z ) — 
{n itm (r, z),t^ m (r, z),V±4>i !m (r, z),v^ m (r, z)}, and c ijm (t) are the complex time-dependent 
amplitudes. The 9 dependence of the basis functions has been explicitly imposed as a 
Fourier basis. The total number of linearly independent basis functions is the number of 
total grid points used in the simulation times the number of independent fields, which is four 
in this case. Now, ipi,m{r,z) can be any linearly independent set of functions and need not 
be the linear eigenfunctions of the system. In fact, the linear eigenfunctions of the equations 
used here are not orthogonal, and are thus not very useful to consider. However, it is quite 
useful to set ipo,m( r , z ) to the fastest growing linear eigenmode because this is the structure 
of interest that is to be projected onto the turbulence. The other VVo,m( r J z ) comprise the 
remainder of the orthonormal basis, and they must be different from the remaining linear 
eigenfunctions in order to complete the orthogonal basis. It isn't necessary for the purpose 
of this study to actually compute these other basis functions, but if one were to compute 
them, one might start with all of the linear eigenmodes and perform a Gram-Schmidt or- 
thogonalization procedure, making sure to start with the fastest growing eigenmode in order 
to preserve it. Using this procedure, Hatch et al.— found that a significant fraction (~ 50%) 
of the energy in a turbulent state of ITG turbulence was contained in the fastest growing 
linear eigenmode at each perpendicular wavenumber. Such a result, however, doesn't require 
knowledge of the other basis functions, and thus they are not computed here. 

Now, to compute this fraction, first define an inner product that is energetically mean- 




(21) 
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ingful and that sets the orthonormality of the basis functions: 



(V>i,m, 1pj,m) = J W^l m * ^j,m^ = (22) 

The weighting w is such that (fturb, fturb) = E turb . Now from Eqs.[2T]and[22], (f turb , f turb ) = 
Eturb = Y.i, m \ c i,™\ 2 and (fturb, m , fturb, m) = E tU rb,m = Y,i\ c i,m\ 2 ■ Then, the amount of en- 
ergy contained in the fastest growing mode (for each m) is given by the square of the 
projection of the mode onto the turbulence: E 0m = \(ipo,m, fturb,m}\ 2 = |co,m| 2 - The ratio 
R m = E 0tm /E tU rb,m is a measure of the fraction of turbulent energy contained in the fastest 
growing linear eigenmode. 

Of course, E turb>m is easily calculated from the turbulent state, but E Q m in the turbulent 
state can only be found with knowledge of the fastest growing eigenf unction. The fastest 
growing eigenfunction, though, can be found easily by running a simulation from a random or 
turbulent state with all of the nonlinearities removed from the model equations as was done 
m Section 1111 After some time, the fastest growing eigenfunctions will come to dominate 
the fluctuation structure. Then, a Fourier decomposition in m space will separate the 
fastest growing eigenfunctions at each m, including the real and imaginary part of the 
eigenfunctions (up to a time dependent complex constant, which is removed by normalizing 
the eigenfunction). These eigenfunctions can then be projected onto the turbulent state 
with the inner product defined in Eq. [22J, giving E 0jm . 

The ratio R m is shown in Fig.[6]for the five simulations. For the most part, the simulations 
other than the n = suppressed one have a small value of the ratio (R m < 0.25) for all m. 
This confirms that the turbulence largely self-organizes without regard to the linear physics. 
The one exception is the Dirichlet simulation for m > 50, which has R m ~ 0.5. Most of 
the energy in this and the other simulations, however, is at low m where the turbulent 
growth rates are largest^, so these larger m eigenmodes don't make a large impact on the 
overall structure of the turbulence. In fact, R m is below 0.1 for m < 40 for the periodic, 
Dirichlet, and Neumann simulations, precisely the area where n = structures dominate the 
energy spectrum^. It is not surprising then that the fastest growing linear eigenfunctions, 
which have little or no flute character, don't make up much of the energy of the signal. 
The sheath simulation shows more linear eigenmode dominance at low m because of the 
relatively large CWM growth rate, but R m is still only about 0.2, which isn't enough for 
the linear eigenmodes to dominate the turbulent structure. Also not surprising is that the 
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fastest growing eigenfunctions make up a significant fraction of the energy in the n = 
suppressed simulation. Where the linear drift wave instability (and the turbulent growth 
rate) is the strongest (at m ~ 50), R m ~ 0.5. The linear physics controls the n = 
suppressed simulation, and the linear eigenmode structure certainly asserts itself in the 
turbulence, but only to a certain degree (50%). Overall, Fig. O shows the relative strength 
of the linear instabilities compared to nonlinear effects for each of the simulations. While the 
nonlinear instability is dominant for all of the simulations other than the n = suppressed 
one, the linear instabilities do still act to varying degree in all of the simulations. 




FIG. 6. The ratio R m of the turbulent energy in the fastest growing linear eigenmode to the total 
energy at each m for all of the simulations. 



VII. CONCLUSION 

The observation of filamentary kn = structures is common in many different kinds 
of experiments and simulations^ - —. Not surprisingly, the presence of these structures is 
usually attributed to linear flute-like instabilities such as flow-driven Kelvin-Helmholtz or 
interchange instabilities rather than to nonlinear instabilities^— . Due to natural limita- 
tions in plasma turbulence experiments, one usually has to resort to indirect evidence to 
gain insight into physical mechanisms of observed turbulence. On the other hand, numerical 
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simulations have all spatial and temporal information available which allows one to per- 
form detailed turbulence analyses, such as the energetics analysis performed in this paper. 
This can make turbulence simulations (provided they are validated on observable data) an 
important tool for uncovering underlying physical mechanisms. 

Through simulation of a particular LAPD experiment, this paper has shown a nonlinear 
instability to be the driving force behind the turbulence. More acurately, this paper has 
extended earlier work— that showed this. However, this extension is important because it 
deals with the matter of axial boundary conditions, and the nonlinear instability depends 
on axial wave dynamics, so the boundary conditions could have greatly affected this. And 
although the various boundary conditions used here do have significant qualitative and/or 
quantitative effects on the linear instabilities of the system, they do not affect the turbulence 
in a significant way. Specifically, the nonlinear instability that preferentially drives k\\ = 
structures remains robust to the change in boundary conditions. Quantitatively, the sheath, 
zero-value, and zero-derivative boundary conditions do cause the linear instability to be 
more competitive with the nonlinear instability, but this effect is not large enough to change 
the qualitative picture. 
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Appendix A: Non-periodic Fourier Decomposition 

It is well known that Fourier reconstructions of signals that contain discontinuities or 
non-periodic boundaries are subject to Gibbs phenomena. A clear indication of this is 
the convergence properties of the Fourier reconstructions. Take a discrete signal with the 
following Fourier decomposition: 



N 




(Al) 



k=-N 
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FIG. 7. D n for a) The simulation with periodic axial boundaries displaying exponential convergence 
and b) the simulation with sheath axial boundaries displaying algebraic convergence. 

where the fk are ordered in the sum by the size of their absolute value with f being the 
largest Fourier coefficient. The Fourier reconstruction of order n < N is then: 



9n(x) = 



(A2) 



k=—n 



There are several types of convergences of the g n , one of which is the LI norm. Defining 
the difference between the original signal and the Fourier reconstruction of order n as D n = 
XL \f( x ) ~ 9n{%)\) one can look at the convergence of D n as a function of n. For periodic 
signals, D n converges exponentially, while it only converges algebraically (power law) for 
non-periodic or discontinuous signals. 

D n is plotted for the cases of the periodic and sheath simulations in Fig. [71 Even though 
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the x-axis label n indicates the mode with the n largest amplitude by construction of 
Eq. \A2\ it also happens to correspond to the axial mode number for all but the last few 
n. In other words, in both simulations, most of the energy is contained in n = modes 
followed by n = 1 modes and so on. Therefore, in reality, the axial Fourier decomposition 
is a useful tool for even the sheath simulation despite the fact that Fourier modes are not 
natural eigenmodes in this case. 



REFERENCES 

1 P. Manneville, Pramana 70, 1009 (2008). 
2 R. E. Waltz, Phys. Rev. Lett. 55, 1098 (1985). 
3 B. D. Scott, Phys. Rev. Lett. 65, 3289 (1990). 
4 B. D. Scott, Phys. Fluids B 4, 2468 (1992). 

5 H. Nordman, V. P. Pavlenko, and J. Weiland, Phys. Fluids B 5, 402 (1993). 
6 K. Itoh, S. I. Itoh, M. Yagi, and A. Fukuyama, Journal of the Physical Society of Japan 
65, 2749 (1996). 

7 E. G. Highcock et al, Zero- Turbulence Manifold in a Toroidal Plasma, Phys. Rev. Lett, 
(submitted). 

8 A. M. Dimits, G. Bateman, and et al., Phys. Plasmas 7, 969 (2000). 

9 D. R. Ernst, P. T. Bonoli, and et al, Phys. Plasmas 11, 2637 (2004). 
10 D. Biskamp and A. Zeiler, Phys. Rev. Lett. 74, 706 (1995). 
11 J. F. Drake, A. Zeiler, and D. Biskamp, Phys. Rev. Lett. 75, 4222 (1995). 
12 A. Zeiler, D. Biskamp, J. F. Drake, and P. N. Guzdar, Phys. Plasmas 3, 2951 (1996). 
13 A. Zeiler, J. F. Drake, and D. Biskamp, Phys. Plasmas 4, 991 (1997). 
14 S. B. Korsholm, P. K. Michelsen, and V. Naulin, Phys. Plasmas 6, 2401 (1999). 
15 B. D. Scott, New J. Physics 4, 52.1 (2002). 
16 B. D. Scott, Plasma Phys. Control. Fusion 45, A385 (2003). 
17 B. D. Scott, Phys. Plasmas 12, 062314 (2005). 
18 B. Friedman et al, Phys. Plasmas 19, 102307 (2012). 
19 W. Gekelman et al, Rev. Sci. Inst. 62, 2875 (1991). 

20 H. L. Berk, D. D. Ryutov, and Y. A. Tsidulko, Phys. Fluids B 3, 1346 (1991). 
21 D. A. Schaffner et al, Phys. Rev. Lett. 109, 135002 (2012). 

22 



S. I. Braginskii, in Reviews of Plasma Physics, edited by M. A. Leontovich (Consultants 

Bureau, New York, ADDRESS, 1965), Vol. 1, pp. 205-311. 

B. D. Dudson et al, Computer Physics Communications 1467 (2009). 

P. Popovich, M. V. Umansky, T. A. Carter, and B. Friedman, Phys. Plasmas 17, 102107 

(2010). 

P. Popovich, M. V. Umansky, T. A. Carter, and B. Friedman, Phys. Plasmas 17, 122312 
(2010). 

M. V. Umansky et al, Phys. Plasmas 18, 055709 (2011). 

B. Friedman, M. V. Umansky, and T. A. Carter, Contrib. Plasma Phys. 52, 412 (2012). 

J. Loizu, P. Ricci, F. Halpern, and S. Jolliet, Phys. Plasmas 13, 122307 (2012). 

D. R. Hatch et al, Phys. Rev. Lett. 106, 115003 (2011). 

B. Rogers and P. Ricci, Phys. Rev. Lett. 104, 225002 (2010). 

K. Kamataki et al, J. Physical Society of Japan 76, 054501 (2007). 

F. Brochard, E. Gravier, and G. Bonhomme, Phys. Plasmas 12, 062104 (2005). 



23 



